"""
Bugs in FermiST
# 1.2.1: Wrong fval in UpperLimit, avoid use this version
"""
from __future__ import print_function
import os
import pprint
import numpy as np
import UnbinnedAnalysis as UbAn
import BinnedAnalysis as BAn
import pyLikelihood as pyLike
from SummedLikelihood import SummedLikelihood
from UpperLimits import UpperLimits
from LikelihoodState import LikelihoodState
from XmlTools import SMKernel

new_fermitools = hasattr(BAn, 'BinnedConfig')
#new_fermitools = False # due to the bug in upperlimit of FermiST 1.2.1, keep back

_small_val = 1.e-6
_one_pls = 1. + _small_val
_one_neg = 1. - _small_val

_irfs_prefix_support_edisp = ['P8R2_', 'CALDB', 'P8R3_']
_iem_prefix = ['gll_iem_', 'gal_2yearp7v6_', 'iem_', 'myiem_']
_iso_prefix = ['iso_', 'myiso_']

_no_covar_algs = ['drmnfb', 'drmngb']

class _MyLikeBase(object):

    def __init__(self, irfs):
        self.irfs = irfs

        self.obs = None
        self.fitobj = None
        self.optobj = None

        self.emin, self.emax = None, None
        self.edisp_flag = None
        self.config = None

        self.centerSrc = None
        self.deleted_snames = []
        self.deleted_srcs = {}

        self.pristines = []

        self._fitmodel = None

    def _checkfitobj(self):
        if self.fitobj is None:
            raise RuntimeError('Please run init fitobj first!')

    def _checkSrcName(self, srcName):
        self._checkfitobj()
        if not isinstance(srcName, str):
            print('[warning] srcName ({}) is not a str'.format(str(srcName)))
        if not srcName in self.fitobj.sourceNames():
            raise IOError('Source %s is not in the model!' % srcName)

    def get_distance(self, srcName):
        """
        Get the angular separation between the input src and centersrc
        """
        self._checkSrcName(srcName)
        if not isinstance(self.centerSrc, str):
            return None
        elif (self.fitobj[self.centerSrc].src.getType() != 'Point') or \
             (self.fitobj[srcName].src.getType() != 'Point'):
            return None
        else:
            return self.fitobj._separation(self.fitobj[self.centerSrc].src, self.fitobj[srcName].src)

    def get_normPar(self, srcName):
        """
        Get the normalization parameter of a source.
        """
        self._checkSrcName(srcName)
        pnames = self.fitobj[srcName].funcs['Spectrum'].paramNames
        if 'Prefactor' in pnames:
            _pname = 'Prefactor'
        elif 'norm' in pnames:
            _pname = 'norm'
        elif 'Integral' in pnames:
            _pname = 'Integral'
        elif 'Value' in pnames:
            _pname = 'Value'
        elif 'Normalization' in pnames:
            _pname = 'Normalization'
        else:
            raise RuntimeError('No invalid normalization for %s!'%srcName)
        return _pname

    def setCenterSrc(self, centerSrc=None):
        if centerSrc is None:
            self.centerSrc = None
        else:
            ctrfile = None
            if centerSrc.startswith('@'):
                ctrfile = centerSrc[1:]
            elif os.path.isfile(centerSrc):
                ctrfile = centerSrc

            if ctrfile is not None:
                srclst = []
                with open(ctrfile, 'r') as fh:
                    for l in fh:
                        if l.startswith('#'):
                            continue
                        lstr = l.strip()
                        try:
                            self._checkSrcName(lstr)
                        except IOError:
                            print('[setCenterSrc] Source %s is not in the model!' % lstr)
                        else:
                            srclst.append(lstr)
                if len(srclst) == 0:
                    pass
                elif len(srclst) == 1:
                    self.centerSrc = srclst[0]
                else:
                    self.centerSrc = srclst
            else:
                try:
                    self._checkSrcName(centerSrc)
                except IOError:
                    print('[setCenterSrc] Source %s is not in the model!' % centerSrc)
                else:
                    self.centerSrc = centerSrc

    def freezeSource(self, srcName, freezeAll=True, disableEdisp=False):
        """
        This is a method to freeze sources. Set freezeAll=True to freeze all the
        parameters of the source. Set disableEdisp=True to disable EDISP for the
        sources of whom all the paramters are fixed.
        """
        self._checkSrcName(srcName)

        if self.fitobj[srcName].src.fixedSpectrum():
            return False

        _freezed = False
        for par in self.fitobj[srcName].funcs['Spectrum'].paramNames:
            if (par == 'Prefactor' or par == 'norm' or
                par == 'Integral' or par == 'Value' or
                par == 'Normalization'):
                normpar = par
            elif self.fitobj[srcName].funcs['Spectrum'].params[par].isFree():
                index = self.fitobj.par_index(srcName, par)
                self.fitobj.freeze(index)
                _freezed = True

        if freezeAll or (not _freezed):
            print('[freezeSrc] All the parameters of %s are freezed' % srcName)
            index = self.fitobj.par_index(srcName, normpar)
            self.fitobj.freeze(index)
        else:
            print('[freezeSrc] The spectrum of %s is freezed' % srcName)

        if self.edisp_flag and disableEdisp:
            try:
                self.fitobj[srcName].src.set_edisp_flag(False)
                print('[freezeSrc] The EDISP of %s is OFF' % srcName)
            except AttributeError:
                pass

        self.init_optobj()
        return True

    def deleteSource(self, srcName):
        """
        Delete a source by setting the scale of normalization to a very
        small number.
        """
        normPar = self.get_normPar(srcName)

        parobj = self.fitobj[srcName].funcs['Spectrum'].params[normPar].parameter
        parobj.setScale(parobj.getScale()*1.e-10)
        self.freezeSource(srcName, freezeAll=True, disableEdisp=True)
        print ('[deleteSource] The %s of %s is scaled down!'%(normPar, srcName))

    def unLoadSource(self, srcName):
        """
        This function removes a source from the model and stores it so that
        you can use it later if you would like.  This is useful if you are
        working on an upper limit and need to get a fit to work before you can
        calculate the upper limit.
        """
        self._checkSrcName(srcName)
        self.deleted_snames.append(srcName)
        self.deleted_srcs[srcName] = self.fitobj.deleteSource(srcName)
        print('[unLoadSource] %s is removed in self.fitobj' % srcName)

    def reLoadSource(self, srcName=None):
        """
        This function puts the last source removed by the unLoadSource function
        back into the model.
        """
        if srcName is None and len(self.deleted_snames) > 0:
            srcName = self.deleted_snames.pop()
        elif srcName in self.deleted_snames:
            self.deleted_snames.remove(srcName)
        else:
           raise RuntimeError("Saved Source does not exist.")

        if self.config is not None:
            self.fitobj.addSource(self.deleted_srcs[srcName], self.config)
        else:
            self.fitobj.addSource(self.deleted_srcs[srcName])
        del self.deleted_srcs[srcName]
        print('[reLoadSource] Reload the source %s.' % srcName)

    def calcUL(self, srcName, ULemin=None, ULemax=None):
        """
        Calculates an upper limit for a source in your model.
        """
        self._checkSrcName(srcName)

        self.ulobj = UpperLimits(self.fitobj)
        if isinstance(ULemin, float) or isinstance(ULemin, int):
            _emin = float(ULemin)
        else:
            if self.emin is None:
                _emin = self.fitobj.energies.min()
            else:
                _emin = self.emin

        if isinstance(ULemax, float) or isinstance(ULemax, int):
            _emax = float(ULemax)
        else:
            if self.emax is None:
                _emax = self.fitobj.energies.max()
            else:
                _emax = self.emax

        print('Calculating 0.95 Upper Limit flux between [%.3f, %.3f]...' % (_emin, _emax))
        print('# IDX  Xval  L-L0  Flux')
        self.ulobj[srcName].compute(emin=_emin, emax=_emax)
        print('The upper limit of %s is %s' % (srcName, self.ulobj[srcName].results))
        return self.ulobj[srcName].results[0].value, _emin, _emax

    def pokeSource(self, srcName, pname=None):
        """
        This function pokes a paramter of a source to a value that is 10%
        of what it was.  This is a useful function to use when you are trying
        to get convergence on a fit.  Many times, a fit won't converge because
        the initial model is too close to the final answer (ie. the minimizer
        does not have enough flexibility to accurately calculate a correlation
        matrix).  In this case, run this funciton on one of the stronger
        sources in your model and redo the fit.  This function is also useful
        to determine how robust your fit is.
        """
        if isinstance(srcName, list):
            srcName = srcName[0]
        if pname is None:
            _pname = self.get_normPar(srcName)
        else:
            self._checkSrcName(srcName)
            if pname in self.fitobj[srcName].funcs['Spectrum'].paramNames:
                _pname = pname
            else:
                raise RuntimeError('Invalid paramName %s!'%pname)

        parobj = self.fitobj[srcName].funcs['Spectrum'].params[_pname].parameter
        old_value = parobj.getValue()
        pbounds = parobj.getBounds()
        if 0.1*old_value < pbounds[0]:
            parobj.setScale(parobj.getScale()*0.1)
        else:
            parobj.setValue(0.1*old_value)
        print ("[pokeSource] Resetting the %s of %s from %.2g to %.2g" %
               (_pname, srcName, old_value, 0.1*old_value))

    def tuneScale(self, srcName=None):
        if srcName is None:
            if self.centerSrc is None:
                raise IOError('Please input source name!')
            else:
                _src = self.centerSrc
        else:
            _src = srcName

        if isinstance(_src, list):
            for __src in _src:
                normpar = self.get_normPar(__src)
                parobj = self.fitobj[__src].funcs['Spectrum'].params[normpar].parameter
                pvalue = parobj.getValue()
                pscale = parobj.getScale()
                newpscale = 10.**(int(np.log10(pvalue*pscale).round()))

                if abs(newpscale/pscale-1.) > _small_val:
                    newpvalue = pvalue*pscale/newpscale
                    parobj.setValue(newpvalue)
                    parobj.setScale(newpscale)
                    print('[tuneScale] [%s__%s] value: %s->%s, scale: %s->%s' % \
                          (__src, normpar, pvalue, newpvalue, pscale, newpscale))
        else:
            normpar = self.get_normPar(_src)
            parobj = self.fitobj[_src].funcs['Spectrum'].params[normpar].parameter
            pvalue = parobj.getValue()
            pscale = parobj.getScale()
            newpscale = 10.**(int(np.log10(pvalue*pscale).round()))

            if abs(newpscale/pscale-1.) > _small_val:
                newpvalue = pvalue*pscale/newpscale
                parobj.setValue(newpvalue)
                parobj.setScale(newpscale)
                print('[tuneScale] [%s__%s] value: %s->%s, scale: %s->%s' % \
                      (_src, normpar, pvalue, newpvalue, pscale, newpscale))

    def pokeBrightest(self, gal_keyword='iem_', iso_keyword='iso_'):
        print('[pokeBrightest] Check the ts value and find the bright source to poke:')
        brightest, TSmax = None, 0.
        for name in self.fitobj.sourceNames():
            if self.fitobj[name].src.fixedSpectrum() or \
               (gal_keyword in name.lower()) or (iso_keyword in name.lower()):
                continue

            ts = self.fitobj.Ts(name)
            print('  [%s] TS=%.3f' % (name, ts))
            if ts > TSmax:
                brightest, TSmax = name, ts

        if brightest is None:
            raise RuntimeError('The fit fails! No reliable information is available.')
        else:
            self.pokeSource(brightest)

    def writeResults(self, outfile='results_pyLike.dat', \
                     ULTSLimit=10., ULemin=None, ULemax=None, NoUpperLimit=False, ApproxTS=True):
        """
        Write the fit results in a dict-like fashion.
        Codes belows are adapted from the one of Y-F, Liang
        ApproxTS=False to refit the TS value of the centerSrc (no effect on other srcs)
        """
        self._checkfitobj()

        dicttot = {'logLikelihood' : self.fitobj.logLike.value()}
        try:
            dicttot['zEDM'] = self.fitobj.optObject.getDistance()
        except AttributeError:
            pass
        try:
            dicttot['zReturnCode'] = self.fitobj.optObject.getRetCode()
        except AttributeError:
            pass
        try:
            dicttot['zFitQuality'] = self.fitobj.optObject.getQuality()
        except AttributeError:
            pass

        if not (self.emin is None or self.emax is None):
            dicttot['zERange'] = (self.emin, self.emax)

        if not self.edisp_flag is None:
            dicttot['zEdispON'] = str(self.edisp_flag)
        else:
            try:
                dicttot['zEdispON'] = str(self.fitobj.logLike.use_edisp())
            except AttributeError:
                try:
                    dicttot['zEdispON'] = str(self.fitobj.logLike.edisp_val()>=0)
                except:
                    dicttot['zEdispON'] = 'N/A'

        lstctrsrc = isinstance(self.centerSrc, list)
        cts_all, cts_free = 0., 0.
        nobs = self.fitobj.total_nobs()
        for src in self.fitobj.sourceNames():
            src_dict = {}

            npred = self.fitobj.NpredValue(src)
            src_dict['NPred'] = npred

            cts_all += npred

            if not self.fitobj[src].src.fixedSpectrum():
                srckey = src
                cts_free += npred

                if (self.centerSrc is not None) and (not ApproxTS) and \
                        ((not lstctrsrc and src == self.centerSrc) or
                                (lstctrsrc and src in self.centerSrc)):
                    print('[INFO] calc accurate TS for {}'.format(src))
                    ts = self.fitobj.Ts(src, reoptimize=True)
                else:
                    ts = self.fitobj.Ts(src)
                src_dict['TS value'] = ts
                line = '  [%s] TS=%.3f, npred=%.1f' % (src, ts, npred)

                if any([src.startswith(iem_pfx) for iem_pfx in _iem_prefix]):
                    normpar_iem = self.get_normPar(src)
                else:
                    normpar_iem = None

                if any([src.startswith(iso_pfx) for iso_pfx in _iso_prefix]):
                    normpar_iso = self.get_normPar(src)
                else:
                    normpar_iso = None

                for pname in self.fitobj[src].funcs['Spectrum'].paramNames:
                    pvalue = self.fitobj[src].funcs['Spectrum'].params[pname].value()
                    if self.fitobj[src].funcs['Spectrum'].params[pname].isFree():
                        perror = self.fitobj[src].funcs['Spectrum'].params[pname].error()
                        src_dict[pname] = '%f +/- %f' % (pvalue, perror)
                    else:
                        src_dict[pname] = '%f' % pvalue

                    if (not normpar_iem is None) and pname == normpar_iem:
                        line += ', %s=%.3f' % (normpar_iem, pvalue)

                    if (not normpar_iso is None) and pname == normpar_iso:
                        line += ', %s=%.3f' % (normpar_iso, pvalue)

                try:
                    if self.emin is None:
                        myemin = 100.
                        print(' [INFO] use the lower bound of 100 MeV to calc flux')
                    else:
                        myemin = self.emin
                    if self.emax is None:
                        myemax = 1.e6
                        print(' [INFO] use the upper bound of 1e6 MeV to calc flux')
                    else:
                        myemax = self.emax
                    flux = self.fitobj.flux(src, emin=myemin, emax=myemax, energyFlux=False)
                    eflux = self.fitobj.flux(src, emin=myemin, emax=myemax, energyFlux=True)
                    if self.fitobj.covar_is_current:
                        fluxErr = self.fitobj.fluxError(src, emin=myemin, emax=myemax, energyFlux=False)
                        efluxErr = self.fitobj.fluxError(src, emin=myemin, emax=myemax, energyFlux=True)
                        src_dict['Flux'] = '%.5e +/- %.5e' % (flux, fluxErr)
                        src_dict['EFlux'] = '%.5e +/- %.5e' % (eflux, efluxErr)

                        if fluxErr > 0.5*flux or ts < 10. or npred < 2.:
                            check_line = '[W]'
                        else:
                            check_value = flux/(fluxErr*np.sqrt(npred))
                            src_dict['CheckValue'] = check_value
                            if check_value <= 0.5:
                                check_line = '[L!]'
                            elif check_value <= 0.9:
                                check_line = '[L?]'
                            elif check_value <= 1.1:
                                check_line = '[G]'
                            else:
                                check_line = '[H!]'
                    else:
                        src_dict['Flux'] = '%.5e' % flux
                        check_line = ''
                except RuntimeError:
                    check_line = ''

                if (not self.centerSrc is None):# and (self.centerSrc in self.fitobj.sourceNames()):
                    dist = self.get_distance(src)
                    if not dist is None:
                        src_dict['distance'] = '%.3f' % dist
                        line += ', distance=%.3f' % dist

                    line += check_line
                    if (lstctrsrc and src in self.centerSrc) or (not lstctrsrc and src == self.centerSrc):
                        ul_flag1 = (not ULTSLimit is None) and ts <= ULTSLimit
                        ul_flag2 = self.fitobj.covar_is_current and flux < 2.*fluxErr
                        ul_flag3 = npred <= max(1., nobs*1.e-4)
                        if (not NoUpperLimit) and (ul_flag1 or ul_flag2 or ul_flag3):
                            flux_ul, _emin, _emax = self.calcUL(src,\
                                                                ULemin=ULemin, ULemax=ULemax)
                            src_dict['Upper Limit'] = flux_ul
                            src_dict['ULERange'] = (_emin, _emax)
                            src_dict['ULFLAG'] = ul_flag1 + 2*ul_flag2 + 4*ul_flag3

                        line += ' [CENTER SOURCE]'
                else:
                    line += check_line
                print(line)
            else:
                srckey = '+'+src
                if npred >= 1.e-4*nobs and npred >= 2.:
                    ts = self.fitobj.Ts(src)
                    src_dict['NPred'] = npred
                    src_dict['TS value'] = ts
                    if npred >= 1.e-2*nobs or ts >= 400.:
                        print('  [%s] TS=%.3f, npred=%.1f [FIXED]' % (src, ts, npred))

            dicttot[srckey] = src_dict

        dicttot['zCounts'] = (cts_all, cts_free, nobs)
        fh = open(outfile, 'w')
        pprint.pprint(dicttot, stream=fh)
        fh.close()
        print ("[writeResults] Counts: Model(Free)/Observed = %.2f(%.2f)/%.1f = 1.(%.4f)%+.6f" %
               (cts_all, cts_free, nobs, cts_free/cts_all, cts_all/nobs-1.))

    def paramsAtLimit(self):
        """
        This function will print out any sources whose parameters
        are close to their limits.  You could use this to find sources
        that are having issues being fit.  This function is useful
        when you're having trouble getting convergence from the New
        Minuit fit routine.
        """
        self._checkfitobj()

        for src in self.fitobj.sourceNames():
            for name in self.fitobj[src].funcs['Spectrum'].paramNames:
                parobj = self.fitobj[src].funcs['Spectrum'].getParam(name)
                if parobj.isFree():
                    pmin, pmax = sorted(parobj.getBounds())
                    pval = parobj.value()
                    perr = parobj.error()

                    lowlimit_flag, upplimit_flag = False, False

                    if (pmin > 0) and (pmax > 100.*pmin):
                        limit_fac = (pmax/pmin)**0.2
                        if pval < pmin * limit_fac:
                            lowlimit_flag = True
                        if pval > pmax / limit_fac:
                            upplimit_flag = True
                    else:
                        limit_fac = (pmax-pmin)*0.001
                        if pval < pmin + limit_fac:
                            lowlimit_flag = True
                        if pval > pmax - limit_fac:
                            upplimit_flag = True

                    if lowlimit_flag:
                        print ('[paramsAtLimit] %s[%s]: value (%s) at its lower limit (%s)!' %
                               (src, name, pval, pmin))
                        continue
                    if upplimit_flag:
                        print ('[paramsAtLimit] %s[%s]: value (%s) at its upper limit (%s)!' %
                               (src, name, pval, pmax))
                        continue

                    if pval - perr <= pmin:
                        print ('[paramsAtLimit] %s[%s]: value(%s)-error(%s)<=min(%s)!' %
                               (src, name, pval, perr, pmin))
                        continue
                    if pval + perr >= pmax:
                        print ('[paramsAtLimit] %s[%s]: value(%s)+error(%s)>=max(%s)!' %
                               (src, name, pval, perr, pmax))

    def writeCountsSpectra(self, outfile='counts_spectra_pyLike.fits'):
        self._checkfitobj()
        self.fitobj.writeCountsSpectra(outfile)
        print('[writeCountsSpectra] Save counts spectra as %s!' % outfile)

    def writeCovar(self, outfile='covariance_pyLike.dat',
                   corr_limit1=0.8, corr_limit2=0.6):
        self._checkfitobj()
        lstctrsrc = isinstance(self.centerSrc, list)
        if self.fitobj.covar_is_current:
            fh = open(outfile, 'w')
            fh.write('-'*20+'free parameters'+'-'*20 + '\n')
            iprm, prm_lst = 1, []
            prmc_index = [] # record the index of the center source parameters
            for prm in self.fitobj.params():
                if prm.parameter.isFree():
                    freePrm = '%s[%s]' % (prm.srcName, prm.parameter.getName())
                    fh.write('%s\t%s\n' % (iprm, freePrm))
                    prm_lst.append(freePrm)
                    if (not self.centerSrc is None) and \
                            ((lstctrsrc and prm.srcName in self.centerSrc) or
                                    (not lstctrsrc and prm.srcName == self.centerSrc)):
                        prmc_index.append(iprm-1)
                    iprm += 1

            covar = np.array(self.fitobj.covariance)
            fh.write('-'*20+ 'eigen values of Hessian mattrix'+ '-'*20 + '\n')
            try:
                eigvals_of_hessian = 1./np.linalg.eigvals(covar)
                np.savetxt(fh, eigvals_of_hessian, delimiter='  ', fmt='% .5e')
            except:
                fh.write('Fail!\n')

            fh.write('-'*20+ 'covariance matrix'+ '-'*20+'\n')
            np.savetxt(fh, covar, delimiter='  ', fmt='% .3e')

            fh.write('-'*20+ 'correlation matrix'+ '-'*20 + '\n')
            cov_diag = covar.diagonal()
            if (cov_diag < 0.).any():
                corr = None
                fh.write('covar.diagonal < 0\n')
                print('[writeCovar] covar.diagonal < 0')
            else:
                sigma = np.sqrt(cov_diag)
                corr = covar / (sigma * sigma.reshape(-1, 1))
                np.savetxt(fh, corr, delimiter='  ', fmt='% .3f')

                fh.write('='*20+ 'correlation check'+ '='*20+'\n')
                triu = np.abs(np.triu(corr, 1))

                fh.write('.'*9+ 'Very strong correlation (|corr| >= %f)'%corr_limit1+ '.'*9 + '\n')
                indics = (triu >= corr_limit1).nonzero()
                for i, j in zip(indics[0], indics[1]):
                    line = '%s <-> %s : %s' % (prm_lst[i], prm_lst[j], corr[i, j])
                    fh.write(line + '\n')
                    print('[writeCovar] %s [VERY STRONG CORR]' % line)

                fh.write('.'*9+ \
                         'Strong correlation (%f <= |corr| < %f)'%(corr_limit2, corr_limit1)+ \
                         '.'*8 + '\n')
                indics = ((triu < corr_limit1)*(triu >= corr_limit2)).nonzero()
                for i, j in zip(indics[0], indics[1]):
                    line = '%s <-> %s : %s' % (prm_lst[i], prm_lst[j], corr[i, j])
                    fh.write(line + '\n')
                    print('[writeCovar] %s [STRONG CORR]' % line)
            fh.close()
            print('[writeCovar] Save convariance check as %s!' % outfile)

            ## make a npz file
            outdata = dict(pnames=prm_lst, cov=covar)
            if corr is not None:
                outdata['corr'] = corr
            outnpz = '.'.join(outfile.split('.')[:-1])+'.npz'
            np.savez_compressed(outnpz, **outdata)
        else:
            print('[writeCovar] No covariance is available!')

    def removeWeak(self, tslimit=0, distlimit=0, RemoveFree=False, RemoveFixed=False):
        """
        This function has two main uses: it will print out details
        on all of the sources in your model and it will remove sources
        according to different requirements.  If you just want to
        print out details, execute it this way:

        <obj>.removeWeak(<my_source>)

        Where <obj> is the quickLike object you're using here and
        <my_source> is the name of your source of interest.  You can
        then remove some of these sources from the model if you like.
        For example, if you want to remove all of the fixed sources
        with TS values less than 1, execute it this way:

        <obj>.removeWeak(<my_source>,tslimit=1,RemoveFixed=True)

        You can mix and match any of the options.  You could remove
        all sources (fixed and free) that are below a TS value of 3
        and are 10 degrees from your source of interest by executing:

        <obj>.removeWeak(<my_source>,tslimit=3,distlimit=10,RemoveFree=True,RemoveFixed=True)
        """
        self._checkfitobj()

        flag = False
        for sname in self.fitobj.sourceNames():
            if self.fitobj[sname].src.fixedSpectrum() and (not RemoveFixed):
                continue

            remove = False
            sourceTS = self.fitobj.Ts(sname)
            distance = self.get_distance(sname) or 0.

            if self.fitobj[sname].src.fixedSpectrum():
                free = False
                if (sourceTS < tslimit) and (distance > distlimit) and RemoveFixed:
                    remove = True
            else:
                free = True
                if (sourceTS < tslimit) and (distance > distlimit) and RemoveFree:
                    remove = True

            if remove:
                print("[removeWeak] Removing %s, TS: %s, Free: %s, Distance: %s" % \
                      (sname, sourceTS, free, distance))
                self.deleteSource(sname)
                flag = True
            else:
                print("[removeWeak] Retaining %s, TS: %s, Free: %s, Distance: %s" % \
                      (sname, sourceTS, free, distance))
        return flag

    def freezeWeak(self, tsLimit0=1., tsLimit1=10., tsLimit2=25., distLimit=0.):
        """
        The sources with ts-value < tsLimit0 will be deleted, the sources with
        tsLimit0 < ts-value < tsLimit1 and distance > distLimit will be totally fixed,
        between tsLimit1 and tsLimit2 spectrum will be fixed.
        tsLimit=None will remove the ts limit.
        Edisp will not be changed!
        """
        self._checkfitobj()

        flag = False
        lstctrsrc = isinstance(self.centerSrc, list)
        for name in self.fitobj.sourceNames():
            if self.fitobj[name].src.fixedSpectrum():
                continue
            if (not self.centerSrc is None) and (
                    (lstctrsrc and name in self.centerSrc) or
                    (not lstctrsrc and name == self.centerSrc)):
                continue

            ts = self.fitobj.Ts(name)
            dist = self.get_distance(name)
            if dist is None:
                dist = 0.1
                strings = '  [%s] TS=%.3f' % (name, ts)
            else:
                strings = '  [%s] TS=%.3f, distance=%.3f' % (name, ts, dist)

            if dist > distLimit:
                if (not tsLimit0 is None) and ts < tsLimit0:
                    self.deleteSource(name)
                    strings += ' [DELETED]'
                    flag = True
                elif (not tsLimit1 is None) and ts < tsLimit1:
                    self.freezeSource(name, True, False)
                    strings += ' [ALL-FROZEN]'
                    flag = True
                elif (not tsLimit2 is None) and ts < tsLimit2:
                    self.freezeSource(name, False, False)
                    strings += ' [INDEX-FROZEN]'
                    flag = True
            print(strings)
        return flag

    def saveXml(self, outfile):
        """
        Save xmlfile, if fitted.
        """
        self._checkfitobj()
        self.fitobj.writeXml(outfile)
        print('[saveXml] Save current model as %s!' % outfile)

    def fit(self, covar=True, modelOut='model_pyLike.xml', pokeSrc=None,
            tuneScale=True, checkloglike=True, verbose=True):
        """
        Run model fit.
        covar=True to calc covariance.
        pokeSrc=None, poke the brightest; pokeSrc=False, do not poke; poke=srcName, poke this source.
        tuneScale=True, change the scale of the centerSrc if needed.
        """
        if verbose:
            _verbosity = 3
        else:
            _verbosity = 0

        self._checkfitobj()

        tol0 = self.fitobj.tol
      # loglike0 = self.fitobj.logLike.value() # buggy

        _opt = self.fitobj.optimizer.lower()
        if _opt == 'drmnfb' or _opt == 'drmngb' or \
           _opt == 'newminuit' or _opt == 'minuit':
            try:
                if self.optobj is None:
                    if _opt in _no_covar_algs:
                        self.fitobj.fit(covar=False, verbosity=_verbosity)
                    else:
                        self.fitobj.fit(covar=covar, verbosity=_verbosity)
                else:
                    if _opt in _no_covar_algs:
                        self.fitobj.fit(covar=False, optObject=self.optobj, verbosity=_verbosity)
                    else:
                        self.fitobj.fit(covar=covar, optObject=self.optobj, verbosity=_verbosity)

                # useless since loglike0 is not available
              # loglike1 = self.fitobj.logLike.value()
              # if checkloglike and loglike0 - loglike1 > 0.5:
              #     raise RuntimeError('logLike[%s] < logLike0[%s]' % \
              #                        (loglike1, loglike0))
            except RuntimeError as err:
                loglike0 = self.fitobj.logLike.value()
                print('[fit] Initial Fit Failed: ' + str(err).strip())
                if verbose:
                    print("[fit] logLikelihood: %s" % loglike0)
                if modelOut is None:
                    if verbose:
                        print('[fit] modelOut is None! Do not save xmlfile!')
                else:
                    self.saveXml(modelOut)

                try:
                    print("[fit] Trying tighter tolerance (DRMtol*0.1)")
                    if pokeSrc is None:
                        self.pokeBrightest()
                    elif pokeSrc is False:
                        pass
                    else:
                        self.pokeSource(pokeSrc)

                    if self.optobj is None:
                        if _opt in _no_covar_algs:
                            self.fitobj.fit(covar=False, verbosity=_verbosity, tol=tol0*0.1)
                        else:
                            self.fitobj.fit(covar=covar, verbosity=_verbosity, tol=tol0*0.1)
                    else:
                        if _opt in _no_covar_algs:
                            self.fitobj.fit(covar=False, optObject=self.optobj,
                                            verbosity=_verbosity, tol=tol0*0.1)
                        else:
                            self.fitobj.fit(covar=covar, optObject=self.optobj,
                                            verbosity=_verbosity, tol=tol0*0.1)

                    loglike1 = self.fitobj.logLike.value()
                    if checkloglike and loglike0 - loglike1 > 0.5:
                        raise RuntimeError('logLike[%s] < logLike0[%s]' % (loglike1, loglike0))
                except RuntimeError as err:
                    loglike1 = self.fitobj.logLike.value()
                    print('[fit] Second Fit Failed: ' + str(err).strip())
                    if verbose:
                        print("[fit] logLikelihood: %s" % loglike1)
                    if loglike1 >= loglike0:
                        loglike0 = loglike1
                        if modelOut is None:
                            if verbose:
                                print('[fit] modelOut is None! Do not save xmlfile!')
                        else:
                            self.saveXml(modelOut)

                    try:
                        print("[fit] Trying looser tolerance (DRMtol*10.)")
                        if pokeSrc is None:
                            self.pokeBrightest()
                        elif pokeSrc is False:
                            pass
                        else:
                            self.pokeSource(pokeSrc)

                        if self.optobj is None:
                            if _opt in _no_covar_algs:
                                self.fitobj.fit(covar=False, verbosity=_verbosity, tol=tol0*10.)
                            else:
                                self.fitobj.fit(covar=covar, verbosity=_verbosity, tol=tol0*10.)
                        else:
                            if _opt in _no_covar_algs:
                                self.fitobj.fit(covar=False, optObject=self.optobj,
                                                verbosity=_verbosity, tol=tol0*10.)
                            else:
                                self.fitobj.fit(covar=covar, optObject=self.optobj,
                                                verbosity=_verbosity, tol=tol0*10.)

                        loglike1 = self.fitobj.logLike.value()
                        if checkloglike and loglike0 - loglike1 > 0.5:
                            raise RuntimeError('logLike[%s] < logLike0[%s]' % \
                                               (loglike1, loglike0))
                    except RuntimeError as err:
                        loglike1 = self.fitobj.logLike.value()
                        print('[fit] Third Fit Failed: ' + str(err).strip())
                        if verbose:
                            print("[fit] logLikelihood: %s" % loglike1)
                        if loglike1 < loglike0:
                            return False
        else:
            raise RuntimeError('Invalid fitobj.optimizer %s!' % self.fitobj.optimizer)

        retCode = self.fitobj.optObject.getRetCode()
        if verbose:
            print('[fit] logLikelihood: %r' % self.fitobj.logLike.value())
            try:
                print('[fit] Fit Distance:  %r' % self.fitobj.optObject.getDistance())
            except AttributeError:
                pass
            try:
                print('[fit] Fit Quality:   %s' % self.fitobj.optObject.getQuality())
            except AttributeError:
                pass
            print('[fit] Return Code:   %s' % retCode)

        if tuneScale and (not self.centerSrc is None):
            self.tuneScale(self.centerSrc)

        if modelOut is None:
            if verbose:
                print('[fit] modelOut is None! Do not save xmlfile!')
        else:
            self.saveXml(modelOut)
            self._fitmodel = modelOut

        if retCode > 0:
            print('[fit] FIT DID NOT CONVERGE!!!')
            if verbose:
                failure = ""
                if retCode & 1:
                    failure += " IsAboveMaxEdm"
                if retCode & 2:
                    failure += " HasCovariance"
                if retCode & 4:
                    failure += " HesseFailed"
                if retCode & 8:
                    failure += " HasMadePosDefCovar"
                if retCode & 16:
                    failure += " HasPosDefCovar"
                if retCode & 32:
                    failure += " HasAccurateCovar"
                if retCode & 64:
                    failure += " HasValidCovariance"
                if retCode & 128:
                    failure += " HasValidParameters"
                if retCode & 256:
                    failure += " IsValid"
                print('[fit]' + failure)
            return False
        else:
            print('[fit] FIT CONVERGED!')
            return True

    def init_optobj(self):
        """
        Inititiate optobj.
        """
        _opt = self.fitobj.optimizer.lower()
        if _opt == 'newminuit':
            self.optobj = pyLike.NewMinuit(self.fitobj.logLike)
        elif _opt == 'minuit':
            self.optobj = pyLike.Minuit(self.fitobj.logLike)
        else:
            self.optobj = None

    def savefitobj(self):
        """
        Save the state of fitobj.
        DO NOT USE WHEN NO FIT IS RUN!
        """
        self._pristine = LikelihoodState(self.fitobj)
        self.pristines.append(self._pristine)

    def restorelastfitobj(self):
        """
        Restore to the latest state of fitobj.
        It has nothing to do with self.pristines.
        """
        self._pristine.restore()
        self.init_optobj()

    def restorefitobj(self):
        """
        Restore to the state of fitobj.
        """
        assert len(self.pristines) > 0
        pristine = self.pristines.pop()
        pristine.restore()
        self.init_optobj()

    def addPointSource(self, name, ra, dec,
                       Prefactor=1.e-13, Index=-2.4, Scale=1e3,
                       Prefactor_free=True, Index_free=True):
        """
        add a point source with a power law.
        """
        self._checkfitobj()

        _pref_scale = 10.**(int(np.log10(Prefactor).round()))
        _pref_value = Prefactor / _pref_scale
        _idx_scale = -1.*np.sign(Index)
        _idx_value = Index / _idx_scale
        _idx_min, _idx_max = sorted(np.array([1., 5.]) * _idx_scale)

        try:
            newps = pyLike.PointSource(0, 0, self.fitobj.observation.observation)
        except AttributeError:
            newps = pyLike.PointSource(0, 0, self.fitobj.binnedData.observation)

        pl = pyLike.SourceFactory_funcFactory().create('PowerLaw')
        pl.setParamValues((_pref_value, _idx_value, Scale))

        prefobj = pl.getParam('Prefactor')
        prefobj.setBounds(1e-5, 1e3)
        prefobj.setScale(_pref_scale)
        if Prefactor_free:
            prefobj.setFree(True)
        else:
            prefobj.setFree(False)
        pl.setParam(prefobj)

        idxobj = pl.getParam('Index')
        idxobj.setBounds(_idx_min, _idx_max)
        idxobj.setScale(_idx_scale)
        if Index_free:
            idxobj.setFree(True)
        else:
            idxobj.setFree(False)
        pl.setParam(idxobj)

        newps.setName(name)
        newps.setSpectrum(pl)
        newps.setDir(ra, dec, True, False)
        if self.config is None:
            self.fitobj.addSource(newps)
        else:
            self.fitobj.addSource(newps, self.config)

    def loglikeProfile(self, srcName, outfile='loglike_profile_pyLike.dat',
                       max_dlike=20., max_delta_sqdlike=0.1, dofit=False,
                       flux_min=None, flux_max=None):
        """
        Fix all bkg srcs, vary the prefactor of src and get loglike profile.
        We do trials assuming the prefactor satisfies Gaussian profile.
        Get the profile of negLoglike!
        """
        self._checkSrcName(srcName)
        if not self.fitobj.covar_is_current:
            raise RuntimeError('Please run fit first! check covar=True.')

        assert max_dlike > 0
        assert max_delta_sqdlike > 0
        _sqrt_max_delta_sqdlike = np.sqrt(max_delta_sqdlike)

        self.savefitobj()
        self.freezeSource(srcName, freezeAll=True, disableEdisp=False)
        if dofit:
            self.fit(covar=False, modelOut=None, checkloglike=False,
                    tuneScale=False, verbose=False)

        normPar = self.get_normPar(srcName)
        parobj = self.fitobj[srcName].funcs['Spectrum'].params[normPar]
        pval0 = parobj.parameter.getValue()
        perr0 = max(parobj.error(), 0.1*pval0)
        pscale = parobj.parameter.getScale()
        pmin, pmax = parobj.parameter.getBounds()

        if self.emin is None:
            myemin = 100.
            print(' [INFO] use the lower bound of 100 MeV to calc flux')
        else:
            myemin = self.emin
        if self.emax is None:
            myemax = 1e6
            print(' [INFO] use the upper bound of 1e6 MeV to calc flux')
        else:
            myemax = self.emax

        flux0 = self.fitobj.flux(srcName, emin=myemin, emax=myemax, energyFlux=False)
        llike0 = self.fitobj.logLike.value()
        prefactor_list, loglike_list = [pval0], [llike0]

      # recursive function to split parameter
        def getsplitval(pval_min, pval_max, sqllike_min, sqllike_max):
            _pval_trial = (pval_min + pval_max)/2.
            parobj.parameter.setValue(_pval_trial)
            if dofit:
                self.fit(covar=False, modelOut=None, checkloglike=False,
                         tuneScale=False, verbose=False)
            _llike_trial = self.fitobj.logLike.value()
            _sqllike_trial = np.sqrt(abs(_llike_trial-llike0))

            print(_pval_trial, llike0-_llike_trial, \
                  abs(sqllike_min-_sqllike_trial), abs(_sqllike_trial-sqllike_max))

            _pval_trial_list, _llike_trial_list = [_pval_trial], [_llike_trial]
            if (abs(sqllike_min-_sqllike_trial) > 1e-2*_sqrt_max_delta_sqdlike) and \
               (abs(_sqllike_trial-sqllike_max) > 1e-2*_sqrt_max_delta_sqdlike):
                if abs(sqllike_min-_sqllike_trial) > max_delta_sqdlike:
                    _pval_trial_ext, _llike_trial_ext = getsplitval(pval_min, _pval_trial,
                                                                    sqllike_min, _sqllike_trial)
                    _pval_trial_list.extend(_pval_trial_ext)
                    _llike_trial_list.extend(_llike_trial_ext)

                if abs(sqllike_max-_sqllike_trial) > max_delta_sqdlike:
                    _pval_trial_ext, _llike_trial_ext = getsplitval(_pval_trial, pval_max,
                                                                    _sqllike_trial, sqllike_max)
                    _pval_trial_list.extend(_pval_trial_ext)
                    _llike_trial_list.extend(_llike_trial_ext)
            return _pval_trial_list, _llike_trial_list

        # negative
        atlimit = False
        itime, fac = 1., 1.
        _prefactor_list0, _loglike_list0 = [pval0], [llike0]
        _prefactor_list1, _loglike_list1 = [], []

        while True:
            pval = pval0 - np.sqrt(2.) * itime * perr0 * max_delta_sqdlike * fac

            if pval <= pmin:
                pval = pmin * _one_pls
                atlimit = True

            parobj.parameter.setValue(pval)
            if dofit:
                self.fit(covar=False, modelOut=None, checkloglike=False,
                         tuneScale=False, verbose=False)
            llike = self.fitobj.logLike.value()
            print(pval, llike0-llike)

            _prefactor_list0.append(pval)
            _loglike_list0.append(llike)

            sqllike1 = np.sqrt(abs(_loglike_list0[-2]-llike0))
            sqllike2 = np.sqrt(abs(_loglike_list0[-1]-llike0))
            dsqllike = abs(sqllike1 - sqllike2)
            if dsqllike > 1e2 * max_dlike:
                raise RuntimeError('Too large dsqllike (%s)! There may be some problems in the fit!'%dsqllike)
            elif dsqllike > max_delta_sqdlike:
                pval1, pval2 = _prefactor_list0[-2], _prefactor_list0[-1]
                pvals_split, llike_split = getsplitval(pval1, pval2, sqllike1, sqllike2)
                _prefactor_list1.extend(pvals_split)
                _loglike_list1.extend(llike_split)
            elif dsqllike < 0.1 * max_delta_sqdlike:
                fac *= 1.9
                itime /= 2.

            if atlimit:
                break
            elif llike0-llike >= max_dlike:
                if not flux_min is None:
                    flux = self.fitobj.flux(srcName, emin=myemin, emax=myemax,
                                            energyFlux=False)
                    if flux < flux_min:
                        break
                    print('[loglikeProfile] dlike (%s) is reached the limit, ' % (llike0-llike) + \
                          'but flux (%s) is still larger than flux_min (%s)!' % (flux, flux_min))
                else:
                    break
            itime += 1.

        prefactor_list.extend(_prefactor_list0[1:])
        loglike_list.extend(_loglike_list0[1:])
        prefactor_list.extend(_prefactor_list1)
        loglike_list.extend(_loglike_list1)

        # positive
        atlimit = False
        itime, fac = 1., 1.
        _prefactor_list0, _loglike_list0 = [pval0], [llike0]
        _prefactor_list1, _loglike_list1 = [], []

        while True:
            pval = pval0 + np.sqrt(2.) * itime * perr0 * max_delta_sqdlike * fac

            if pval >= pmax:
                pval = pmax * _one_neg
                atlimit = True

            parobj.parameter.setValue(pval)
            if dofit:
                self.fit(covar=False, modelOut=None, checkloglike=False,
                         tuneScale=False, verbose=False)
            llike = self.fitobj.logLike.value()
            print(pval, llike0-llike)

            _prefactor_list0.append(pval)
            _loglike_list0.append(llike)

            sqllike1 = np.sqrt(abs(_loglike_list0[-2]-llike0))
            sqllike2 = np.sqrt(abs(_loglike_list0[-1]-llike0))
            dsqllike = abs(sqllike1 - sqllike2)
            if dsqllike > 1e2 * max_dlike:
                raise RuntimeError('Too large dsqllike (%s)! There may be some problems in the fit!'%dsqllike)
            elif dsqllike > max_delta_sqdlike:
                pval1, pval2 = _prefactor_list0[-2], _prefactor_list0[-1]
                pvals_split, llike_split = getsplitval(pval1, pval2, sqllike1, sqllike2)
                _prefactor_list1.extend(pvals_split)
                _loglike_list1.extend(llike_split)
            elif dsqllike < 0.1 * max_delta_sqdlike:
                fac *= 1.9
                itime /= 2.

            if atlimit:
                break
            elif llike0-llike >= max_dlike:
                if not flux_max is None:
                    flux = self.fitobj.flux(srcName, emin=myemin, emax=myemax,
                                            energyFlux=False)
                    if flux > flux_max:
                        break
                    print('[loglikeProfile] dlike (%s) is reached the limit, ' % (llike0-llike) + \
                          'but flux (%s) is still smaller than flux_max (%s)!' % (flux, flux_max))
                else:
                    break
            itime += 1.

        prefactor_list.extend(_prefactor_list0[1:])
        loglike_list.extend(_loglike_list0[1:])
        prefactor_list.extend(_prefactor_list1)
        loglike_list.extend(_loglike_list1)

      # make a numpy array and save
        prefactor_list = np.atleast_1d(prefactor_list)
        loglike_list = np.atleast_1d(loglike_list)
        flux_list = flux0/pval0*prefactor_list
        dloglike_list = loglike_list - llike0

        argsort = np.argsort(prefactor_list)
        outdata = np.c_[prefactor_list*pscale, flux_list, dloglike_list][argsort, :]
        np.savetxt(outfile, outdata, delimiter='\t', header='Prefactor\tFlux\tloglike%+.6f'%llike0)

      # set to the previous value
      # parobj.parameter.setValue(pval0)
        self.restorefitobj()


class MyBinnedLike(_MyLikeBase):
    """
    Binned Analysis pyLikelihood script.
    """
    def __init__(self, srcmaps, expcube, ltcube, srcmodel, irfs, edisp=True, edisp_bins=-2):
        _MyLikeBase.__init__(self, irfs)
        self.obs = BAn.BinnedObs(srcMaps=srcmaps, expCube=ltcube,
                                 binnedExpMap=expcube, irfs=self.irfs)
        self.srcmaps, self.srcmodel = srcmaps, srcmodel
        self.ltcube, self.expcube = ltcube, expcube

        if any([irfs.startswith(irfp) for irfp in _irfs_prefix_support_edisp]):
            if edisp is True:
                self.edisp_flag = True
                if self.edisp_flag and new_fermitools:
                    self.config = BAn.BinnedConfig(edisp_bins=edisp_bins)
                    print('[__init__] edisp_bins=%s'%edisp_bins)
            elif edisp is False:
                self.edisp_flag = False
            else:
                self.edisp_flag = None
        else:
            print('[edispSwitch] edisp unsupported irfs??')
            self.edisp_flag = None

    def init_fitobj(self, opt='MINUIT', srcModel=None, tol=1.e-5):
        """
        Initializes a New optimizer.
        """
        if not srcModel is None:
            _scmdl = srcModel
        elif not self._fitmodel is None:
            _scmdl = self._fitmodel
        else:
            _scmdl = self.srcmodel

        first_init_flag = self.fitobj is None

        _opt = opt.lower()
        if new_fermitools:
            if _opt == 'drmnfb':
                self.fitobj = BAn.BinnedAnalysis(self.obs, _scmdl, optimizer="DRMNFB", config=self.config)
            elif _opt == 'drmngb':
                self.fitobj = BAn.BinnedAnalysis(self.obs, _scmdl, optimizer="DRMNGB", config=self.config)
            elif _opt == 'minuit':
                self.fitobj = BAn.BinnedAnalysis(self.obs, _scmdl, optimizer='MINUIT', config=self.config)
            elif _opt == 'newminuit':
                self.fitobj = BAn.BinnedAnalysis(self.obs, _scmdl, optimizer='NewMinuit', config=self.config)
        else:
            if _opt == 'drmnfb':
                self.fitobj = BAn.BinnedAnalysis(self.obs, _scmdl, optimizer="DRMNFB")
            elif _opt == 'drmngb':
                self.fitobj = BAn.BinnedAnalysis(self.obs, _scmdl, optimizer="DRMNGB")
            elif _opt == 'minuit':
                self.fitobj = BAn.BinnedAnalysis(self.obs, _scmdl, optimizer='MINUIT')
            elif _opt == 'newminuit':
                self.fitobj = BAn.BinnedAnalysis(self.obs, _scmdl, optimizer='NewMinuit')

        if self.edisp_flag is not None:
            self.fitobj.logLike.set_edisp_flag(self.edisp_flag)

        self.init_optobj()

        print('[init_fitobj] Something about fitobj:')
        print(self.fitobj)
        try:
            print('[init_fitobj] current edispflag=%s!' % (self.fitobj.logLike.use_edisp()))
        except AttributeError:
            print('[init_fitobj] current edispflag=%s!' % (self.fitobj.logLike.edisp_val()>=0))

        if self.edisp_flag:
            edisp_srcs = []
            for src in list(self.fitobj.model.srcNames):
                if self.fitobj[src].src.use_edisp():
                    edisp_srcs.append(src)
            n_edisp_src = len(edisp_srcs)
            if n_edisp_src > 0:
                print('[init_fitobj] %i source(s) with edisp on: %s' % \
                      (n_edisp_src, ', '.join(edisp_srcs)))
            else:
                print('[init_fitobj] no source with edisp on')

        self.fitobj.tol = float(tol)
        if not first_init_flag:
            self.setERange()
            self.savefitobj()

    def writeModelMap(self, outfile, outtype='CCUBE'):
        self.fitobj.writeModelMap(outfile, outtype=outtype)
        print('[writeModelMap] Save model map as %s!' % outfile)

    def _setEBounds(self, bound_min, bound_max):
        """
        Set Energy range with the number of bounds.
        For example, if you want to fit the first energy bin, you should use self._setEBounds(0, 1)
        """
        self._checkfitobj()

        _bound_min, _bound_max = int(bound_min), int(bound_max)
        self.fitobj.selectEbounds(_bound_min, _bound_max)
        try:
            self.emin, self.emax = self.fitobj.emin, self.fitobj.emax
        except AttributeError:
            self.emin, self.emax = self.fitobj.energies[_bound_min], self.fitobj.energies[_bound_max]
        print('[_setEBounds] ENERGY RANGE is [%s:%s]' % (self.emin, self.emax))

    def setERange(self, Emin=None, Emax=None):
        """
        Use Max(Emin, Emin_default) as minimum energy and Min(Emxa, Emax_default)
        as maximum energy in the fit.
        """
        self._checkfitobj()

        if Emin is None:
            if self.emin is None:
                _emin = self.fitobj.energies.min()
            else:
                _emin = self.emin
            emin_str = '!%s'%_emin
        else:
            _emin = self.fitobj.energies[self.fitobj.energies >= (Emin*_one_neg)].min()
            emin_str = str(_emin)

        if Emax is None:
            if self.emax is None:
                _emax = self.fitobj.energies.max()
            else:
                _emax = self.emax
            emax_str = '!%s'%_emax
        else:
            _emax = self.fitobj.energies[self.fitobj.energies <= (Emax*_one_pls)].max()
            emax_str = str(_emax)

        assert _emin < _emax
        self.fitobj.setEnergyRange(_emin, _emax)
        self.emin, self.emax = _emin, _emax
        print('[setERange] ENERGY RANGE is [%s:%s] MeV' % (emin_str, emax_str))


class MyUnbinnedLike(_MyLikeBase):
    """
    Unbinned Analysis pyLikelihood script.
    """
    def __init__(self, gtifile, expmap, ltcube, scfile, srcmodel, irfs):
        _MyLikeBase.__init__(self, irfs)
        self.obs = UbAn.UnbinnedObs(eventFile=gtifile, scFile=scfile,
                                    expMap=expmap, expCube=ltcube,
                                    irfs=self.irfs)
        self.gtifile, self.scfile = gtifile, scfile
        self.exmap, self.ltcube = expmap, ltcube
        self.srcmodel = srcmodel
        self.edisp_flag = False

    def init_fitobj(self, opt='MINUIT', srcModel=None, tol=1.e-5):
        """
        Initializes a New optimizer.
        """
        if not srcModel is None:
            _scmdl = srcModel
        elif not self._fitmodel is None:
            _scmdl = self._fitmodel
        else:
            _scmdl = self.srcmodel

        first_init_flag = self.fitobj is None

        _opt = opt.lower()
        if _opt == 'drmnfb':
            self.fitobj = UbAn.UnbinnedAnalysis(self.obs, _scmdl, optimizer="DRMNFB")
        elif _opt == 'drmngb':
            self.fitobj = UbAn.UnbinnedAnalysis(self.obs, _scmdl, optimizer="DRMNGB")
        elif _opt == 'minuit':
            self.fitobj = UbAn.UnbinnedAnalysis(self.obs, _scmdl, optimizer='MINUIT')
        elif _opt == 'newminuit':
            self.fitobj = UbAn.UnbinnedAnalysis(self.obs, _scmdl, optimizer='NewMinuit')

        self.init_optobj()

        print('[init_fitobj] Something about fitobj:')
        print(self.fitobj)

        self.fitobj.tol = float(tol)
        if not first_init_flag:
            self.setERange()
            self.savefitobj()

    def setERange(self, Emin=None, Emax=None):
        self.emin, self.emax = self.fitobj.energies.min(), self.fitobj.energies.max()
        print('[setERange] DO NOT SUPPORT UNBINNED LIKELIHOOD! ENERGY RANGE is [%f, %f] MeV'%(self.emin, self.emax))
        if not (Emin is None and Emax is None):
            print('[setERange] Useless parameters: Emin=%f [MeV], Emax=%f [MeV]' % (Emin, Emax))


class MySumLike(_MyLikeBase):

    def __init__(self, irfs):
        _MyLikeBase.__init__(self, irfs)
        self.likeobj_components = {}

    @staticmethod
    def _copyParamValue(infile, example, outfile):
        """
        copy the value of parameters in example to infile, and save to outfile
        """
        x0 = SMKernel.SourceModel(infile, raiseErr=True)
        x1 = SMKernel.SourceModel(example, raiseErr=True)
        for sname, src1 in x1.srcList.items():
            if sname in x0.names and (x0[sname].spectrum.type == src1.spectrum.type):
                src0 = x0[sname]
                for pname, par in src1.spectrum.parameters.items():
                    src0.spectrum.parameters[pname].value = par.value
                    if hasattr(par, 'error'):
                        src0.spectrum.parameters[pname]['error'] = par.error
        x0.writeTo(outfile)

    def addComponent(self, likeobj, name=None):
        """
        Initialize the fitobj of the likeobj and add to sumLike class.
        In order to avoid the failure of fit, use MINUIT as optimizer!
        May contain bug when add components using 'for'?
        If name=None, likeobj will be named as comp1, comp2, comp3, ...
        """
        name_list = self.likeobj_components.keys()
        if isinstance(name, str):
            if name in name_list:
                raise IOError('There is same name object in the memory!')
            _name = name
        else:
            icomp = 0
            while ('comp%i' % icomp) in name_list:
                icomp += 1
            _name = 'comp%i'%icomp

        if not (isinstance(likeobj, MyUnbinnedLike) or
                isinstance(likeobj, MyBinnedLike)):
            raise IOError('Invalid likeobj')

        if likeobj in self.likeobj_components.values():
            raise RuntimeError('There is the same one in the memory!')

        if self.irfs != likeobj.irfs:
            raise RuntimeError('IRFs is not the same!')

        self.likeobj_components[_name] = likeobj

    def saveXml(self, outfile):
        self._checkfitobj()
        self.fitobj.writeXml(outfile)
        print('[saveXml] Save current model as %s!' % outfile)

        _outpfx = '.'.join(outfile.split('.')[:-1])
        for name, likeobj in self.likeobj_components.items():
            _outfile = '%s_%s.xml' % (_outpfx, name)
            self._copyParamValue(likeobj.srcmodel, outfile, _outfile)
            print('[saveXml] Copy the model %s to %s!' % (outfile, _outfile))

    def init_fitobj(self, opt=None, srcModel=None, tol=1.e-5):
        """
        It seems that summed-likelihood only support MINUIT algorithm
        In the method, if srcModel is not None, we copy the parameter of
        srcmodel to the new one in the same folder of input model of likeobjs
        (otherwise the filefunction in isotropic component will be wrong!).
        """
        if not srcModel is None:
            _srcmdl = srcModel
            _srcpfx = '.'.join(_srcmdl.split('.')[:-1])
        elif not self._fitmodel is None:
            _srcmdl = self._fitmodel
            _srcpfx = '.'.join(_srcmdl.split('.')[:-1])
        else:
            _srcmdl = None

        first_init_flag = self.fitobj is None

        self.fitobj = SummedLikelihood()
        self.optobj = None

        for name, likeobj in self.likeobj_components.items():
            if _srcmdl is None:
                srcmdl = None
            else:
                srcmdl = '%s_%s.xml' % (_srcpfx, name)
                if not os.path.isfile(srcmdl):
                    self._copyParamValue(likeobj.srcmodel, _srcmdl, srcmdl)

            if opt is None:
                if likeobj.fitobj is None:
                    likeobj.init_fitobj(opt='MINUIT', srcModel=srcmdl, tol=tol)
            else:
                likeobj.init_fitobj(opt=opt, srcModel=srcmdl, tol=tol)

            self.fitobj.addComponent(likeobj.fitobj)

        if not first_init_flag:
            self.setERange()

    def __iadd__(self, likeobj):
        self.addComponent(likeobj)
        return self

    def _checkfitobj(self):
        if self.fitobj is None:
            raise RuntimeError('Please run init fitobj first!')
            if len(self.likeobj_components) < 1:
                raise RuntimeError('No model is added!')

    def setERange(self, Emin=None, Emax=None):
        print('[setERange] DO NOT SUPPORT SUMMED LIKELIHOOD! PLEASE USE Shell FLAGS!')
        if not (Emin is None and Emax is None):
            print('[setERange] Useless parameters: Emin=%f [MeV], Emax=%f [MeV]' % (Emin, Emax))

    def writeCountsSpectra(self, outfile='counts_spectra_pyLike.fits'):
        print('[writeCountsSpectra] DO NOT SUPPORT SUMMED LIKELIHOOD!')
        if not outfile is None:
            print('[writeCountsSpectra] Useless parameters: outfile=%s' % outfile)

    def init_optobj(self):
        self.optobj = None
        print('[init_optobj] Useless in SUMMED LIKELIHOOD!')


def _fit_sequence(likeobj, centerSrc=None, ApproxTS=True, emin=None, emax=None, saveFolder='./', # basic settings
                  twopasses=True, fitalg1='DRMNFB', fitalg2='MINUIT', tol1=1e-3, tol2=1.e-5,     # fit settings
                  modelout1='model_pyLike_pass1.xml', modelout2='model_pyLike_pass2.xml', # output files
                  NoUpperLimit=False, ULTSLimit=10., ULemin=None, ULemax=None,            # upper limit parameters
                  freezeWeak=False, tuneScale=True, TSLimit0=1., TSLimit1=10., TSLimit2=25., distLimit=0.,# autofreeze parameters
                  likeprofile=False, profile_out='loglike_profile_pyLike.dat',
                  profile_maxdlike=20., profile_maxdlikebin=0.1,  profile_fit=False,      # profile likelihood
                  profile_fluxmin=None, profile_fluxmax=None,
                  results_data='results_pyLike.dat', counts_spec='counts_spectra_pyLike.fits',
                  covar_data='covariance_pyLike.dat' # the name of output data file
                 ):
    """
    Support MyUnbinnedLike, MyBinnedLike, MySumLike now!
    """
    assert isinstance(likeobj, _MyLikeBase)

    if not saveFolder:
        _saveFolder = './'
    elif saveFolder.endswith('/'):
        _saveFolder = saveFolder
    else:
        _saveFolder = saveFolder + '/'

    assert os.path.exists(_saveFolder)

    if modelout1 is None:
        _modelout1 = _saveFolder + 'model_pyLike_pass1.xml'
    elif '/' in modelout1:
        _modelout1 = modelout1
    else:
        _modelout1 = _saveFolder + modelout1

    if modelout2 is None:
        _modelout2 = _saveFolder + 'model_pyLike_pass2.xml'
    elif '/' in modelout2:
        _modelout2 = modelout2
    else:
        _modelout2 = _saveFolder + modelout2

    if profile_out is None:
        _profile_out = _saveFolder + 'loglike_profile_pyLike.dat'
    elif '/' in profile_out:
        _profile_out = profile_out
    else:
        _profile_out = _saveFolder + profile_out

    if results_data is None:
        _resultdat = _saveFolder + 'results_pyLike.dat'
    elif '/' in results_data:
        _resultdat = results_data
    else:
        _resultdat = _saveFolder + results_data

    if counts_spec is None:
        _countsSpec = None
    elif '/' in counts_spec:
        _countsSpec = counts_spec
    else:
        _countsSpec = _saveFolder + counts_spec

    if covar_data is None:
        _covardat = _saveFolder + 'covariance_pyLike.dat'
    elif '/' in covar_data:
        _covardat = covar_data
    else:
        _covardat = _saveFolder + covar_data

    print('likeobj: %s' % likeobj.__class__.__name__)

    if twopasses:
        print('-'*50, 'PASS1', '-'*50)
        likeobj.init_fitobj(opt=fitalg1, tol=tol1)
        likeobj.setCenterSrc(centerSrc)
        likeobj.setERange(emin, emax)

        goodfit = likeobj.fit(False, _modelout1, pokeSrc=likeobj.centerSrc, tuneScale=tuneScale)

        if freezeWeak and goodfit:
            hasfreezed = likeobj.freezeWeak(tsLimit0=TSLimit0, tsLimit1=TSLimit1,
                                            tsLimit2=TSLimit2, distLimit=distLimit)
            if hasfreezed:
                print('-'*50, 'PASS1 REFIT', '-'*50)
                likeobj.fit(False, _modelout1, pokeSrc=likeobj.centerSrc, tuneScale=tuneScale)
        elif (not goodfit) and (likeobj._fitmodel is None):
            likeobj.init_fitobj(opt='MINUIT', tol=tol1)
            likeobj.fit(False, _modelout1, pokeSrc=likeobj.centerSrc, tuneScale=tuneScale)

        print('-'*50, 'PASS2', '-'*50)
        likeobj.init_fitobj(opt=fitalg2, tol=tol2)
    else:
        print('-'*50, 'PASS2', '-'*50)
        likeobj.init_fitobj(opt=fitalg2, tol=tol2)
        likeobj.setCenterSrc(centerSrc)
        likeobj.setERange(emin, emax)

    likeobj.fit(True, _modelout2, pokeSrc=likeobj.centerSrc, tuneScale=False)
    likeobj.paramsAtLimit()

    likeobj.writeResults(_resultdat, ULTSLimit=ULTSLimit, ULemin=ULemin, ULemax=ULemax,
                         NoUpperLimit=NoUpperLimit, ApproxTS=ApproxTS)
    likeobj.writeCovar(_covardat)
    if _countsSpec is not None:
        try:
            likeobj.writeCountsSpectra(_countsSpec)
        except:
            pass

    if likeprofile and (not likeobj.centerSrc is None):
        print('-'*43, 'PROFILE  LIKELIHOOD', '-'*43)
        likeobj.loglikeProfile(srcName=likeobj.centerSrc, outfile=_profile_out, dofit=profile_fit,
                               max_dlike=profile_maxdlike, max_delta_sqdlike=profile_maxdlikebin,
                               flux_min=profile_fluxmin, flux_max=profile_fluxmax)
    print('All done!')

def binLikeCLI(srcMaps, expCube, binnedExpMap, srcModel, irfs, edispon=True, edisp_bins=-2, **kwds):
    print('Making likeobj ...')
    likeobj = MyBinnedLike(srcmaps=srcMaps, expcube=binnedExpMap, ltcube=expCube,
                           srcmodel=srcModel, irfs=irfs, edisp=edispon, edisp_bins=edisp_bins)
    _fit_sequence(likeobj, **kwds)

def unbinLikeCLI(eventFile, scFile, expMap, expCube, srcModel, irfs, **kwds):
    print('Making likeobj ...')
    likeobj = MyUnbinnedLike(gtifile=eventFile, expmap=expMap, ltcube=expCube,
                             scfile=scFile, srcmodel=srcModel, irfs=irfs)
    _fit_sequence(likeobj, **kwds)

def sumLikeFunc(*likeobjs, **kwds):
    print('Making likeobj ...')
    irfs = likeobjs[0].irfs
    sumlikeobj = MySumLike(irfs)
    for likeobj in likeobjs:
        sumlikeobj.addComponent(likeobj)
    _fit_sequence(sumlikeobj, **kwds)
